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This  grant  supported  research  by  two  UC  Berkeley  graduate  students  on 
computational  fluid  dynamics  and  materials  science.  Both  are  developing  inno¬ 
vative  numerical  methods  for  flows  around  complex  boundaries  such  as  occur  in 
solidification  from  the  melt.  Research  summaries  for  each  student  follow. 

Ricardo  Cortez  (currently  an  NSF  postdoctoral  fellow  at  Courant  Insti¬ 
tute)  studied  magnetization  methods  for  fluid  flow. 

A  new  formulation  of  3- dimensional  incompressible  fluid  flow  has  been  pre¬ 
sented  recently  with  the  introduction  of  a  new  variable  sometimes  referred  to 
as  ^‘velocity  magnetization”  or  simply  “magnetization.”  A  discretization  of  the 
resulting  equations  produces  a  Hamiltonian  system  which,  in  turn,  leads  to  a 
new  family  of  Lagrangian  numerical  methods  for  the  study  of  incompressible 
flow.  These  numerical  methods  are  essentially  independent  of  spatial  dimension 
and  have  the  property  that  they  preserve  discrete  invariants  associated  with 
the  Hamiltonian.  However  there  are  many  issues  related  to  this  new  formula¬ 
tion  that  require  further  investigation  in  order  to  arrive  at  efficient  numerical 
algorithms  that  may  be  used  in  a  variety  of  situations.  These  issues  include: 

1.  The  design  of  rohxtsi  hybrid  methods  using  voriicity  and  magneiizaiion 
variables.  The  magnetization  variables  can  be  used  in  conjunction  with 
vortex  methods  in  order  to  exploit  the  advantages  that  each  representation 
offers.  Accurate  and  relatively  inexpensive  procedures  to  switch  from  one 
computational  variable  to  the  other  are  required. 

2.  Procedures  to  ensure  that  the  magnetization  remains  optimal  or  nearly 
optimal.  For  a  given  initial  configuration,  magnetization  naturally  tends  to 
evolve  into  allowable  configurations  that  depart  from  the  optimal  one.  In 
three  dimensions,  finding  the  latter  is  a  problem  related  to  finding  minimal 
surfaces  attached  to  a  fixed  frame.  The  accuracy  of  the  numerical  method 
is  related  to  the  distribution  of  magnetization,  making  the  problem  of 
finding  the  optimal  configuration  a  crucial  part  of  the  method. 
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3.  Applicability  of  these  methods  near  boundaries.  The  formulation  of  the 
Navier-Stokes  equations  in  terms  of  magnetization  has  the  potential  of 
suggesting  new  ways  of  dealing  with  creation  of  vorticity  inside  boundary 
layers. 

4.  Desingularization  of  the  equations  of  motion  and  appropriate  fast  summa¬ 

tion  techniqxies.  The  differential  equations  that  govern  the  evolution  of 
the  particle  positions  and  the  magnetization  they  carry  have  high-order 
singularities,  where  n  is  the  dimension  of  the  problem.  The  or¬ 

der  of  the  singularities  has  an  adverse  effect  on  the  speed  up  attained  by 
some  fast  summation  techniques;  therefore,  an  effort  must  be  made  to 
reduce  the  order  of  the  singularities  in  the  equations  and  to  develop  fast 
summation  procedures  that  are  well  suited  for  these  equations. 

5.  Mathematical  theory.  The  relationship  between  parameters  of  the  method 
and  other  necessary  conditions  for  the  convergence  of  the  discrete  variables 
to  their  continuous  counterparts  are  only  conjectured.  There  is  a  need  for 
theory  that  provides  convergence  and  stability  proofs  and  error  bounds 
based  on  discretization  parameters. 

During  the  academic  years  1993”94  and  1994-95,  this  AASERT  grant  sup¬ 
ported  intensive  work  on  some  of  the  above  issues,  continuing  the  work  on  vortex 
methods  which  was  done  under  the  support  of  AFOSR  Grant  No.  FDF49620- 
93-1-0053  during  the  last  few  years.  A  summary  of  Dr.  Cortez’s  results  follows. 

In  two  dimensions,  the  relationship  between  discrete  magnetization  variables 
and  vortices  has  been  established  by  a  rigorous  interpretation  of  magnetization 
as  vortex  dipoles  with  a  prescribed  dipole  moment.  For  a  given  discretization, 
each  magnetization  vector  can  be  replaced  by  a  pair  of  vortices  of  equal  but 
opposite  strength  located  some  small  distance  apart.  Then,  vortices  sufficiently 
“near”  others  can  be  combined  into  a  single  vortex  with  strength  equal  to  the 
sum  of  individual  strengths.  This  collapse  of  vortices  keeps  the  computation 
feasible.  Error  analysis  provides  a  criterion  for  how  near  vortices  should  be  in 
order  to  be  combined.  The  distance  between  the  vortices  is  a  parameter  which 
is  chosen  so  that  the  error  in  this  procedure  is  of  a  prescribed  size. 


Example  of  combination  of  vortices 

The  reverse  procedure  has  also  been  accomplished.  Given  a  set  of  vortices 
of  net  vorticity  zero,  one  can  extract  subsets  to  be  represented  by  a  string 
of  vortex  pairs,  which  are  then  approximated  by  magnetization  vectors.  The 
errors  in  both  procedures  are  shown  to  be  0{h^),  where  h  is  the  dipole  distance. 
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The  ability  to  switch  from  magnetization  to  vorticity  and  back  allows  numerical 
experiments  that  compare  particle  positions  obtained  by  magnetization  methods 
with  those  obtained  by  well  understood  vortex  methods,  and  provides  a  very 
important  first  step  toward  the  making  of  hybrid  methods.  The  numerical 
experiments  show  the  order  of  convergence. 

In  many  physical  situations,  a  magnetization  configuration,  which  initially 
represents  a  given  set  of  vortices  to  a  prescribed  accuracy,  evolves  naturally  in  a 
way  that  endures  accuracy  loss  to  the  point  where  the  magnetization  no  longer 
approximates  the  original  vortices  evolved  in  time.  The  velocity  field  due  to 
a  string  of  magnetization  vectors  connecting  two  vortices  hcis  been  interpreted 
as  an  approximation  to  the  line  integral  that  represents  the  velocity  field  due 
to  the  vortices.  The  loss  of  accuracy  is  now  understood  as  the  deterioration 
of  the  approximate  line  integral  due  to  the  stretching  of  the  curve  along  which 
we  integrate.  This  understanding  has  led  to  various  procedures  that  correct 
the  problem  and  ensure  that  the  magnetization  remains  optimal  by  maintaining 
the  initial  accuracy  over  time.  Some  of  these  procedures  periodically  discard 
the  integration  curve  and  replace  it  with  a  more  convenient  one,  computing  a 
new  magnetization  field  along  the  new  curve.  Others  use  the  original  curve  but 
refine  the  intervals  used  in  the  quadrature  so  as  to  maintain  the  initial  accuracy. 
The  latter  case  is  less  efficient  since  the  number  of  particles  may  grow  without 
bound,  but  it  has  a  natural  extension  to  S-dimensional  problems.  There,  a 
magnetization  vector  is  approximated  by  a  vortex  loop  on  the  plane  normal  to 
the  vector;  these  loops  can  be  split  into  smaller  ones  to  maintain  the  accuracy 
of  the  discretization.  The  first  draft  of  the  work  above  has  been  written  as  part 
of  Mr.  Cortez’s  dissertation. 

Magnetization  variables  are  particularly  well  suited  for  problems  in  which 
vortex  dipoles  are  used.  One  such  problem  is  the  motion  of  an  elastic  membrane 
immersed  in  a  fluid,  where  the  effects  of  elastic  forces  acting  on  the  fluid  can 
be  introduced  as  the  rate  of  change  of  magnetization.  The  forces  that  the 
membrane  imparts  on  the  fluid  over  a  small  time  interval  represent  impulse, 
and  their  effect  can  be  introduced  via  vortex  dipoles.  In  other  words,  the  forces 
represent  an  evolving  magnetization  field  on  the  membrane.  This  application 
has  been  implemented  in  two  dimensions  for  smooth  membranes  and  inviscid, 
incompressible  flow. 

is  thought  of  as  a  starting  point  for  the  understanding  of  magnetiza¬ 
tion  and  for  providing  extensions  to  R^.  In  three  dimensions,  magnetization 
variables  have  the  advantage  over  vortex  methods  of  automatically  yielding 
a  divergence-free  vorticity  field.  This  is  a  source  of  problems  with  3-D  vor¬ 
tex  methods.  Problems  with  3-D  magnetization  methods  include  the  arbitrary 
shape  of  the  loops  when  two  neighboring  loops  must  be  combined,  and  the  need 
to  find  minimal  surfaces  when  a  new  magnetization  representation  of  large  vor¬ 
tex  loops  is  needed.  Computationally  efficient  extensions  of  the  2-dimensional 
results  will  provide  a  promising  arena  for  the  modeling  of  3-dimensional  fluid 
flow. 
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Hans  Johansen  is  a  doctoral  candidate  in  the  Department  of  Mechanical 
Engineering  at  the  University  of  California,  Berkeley.  His  research  focuses  on 
finite  difference  methods  for  simulation  of  unsteady  fluid  dynamics,  especially  in 
bulk  Czochralski  crystal  growth  of  multi-component  semiconductor  substrates. 
One  of  the  major  factors  that  limits  the  quality  of  such  crystals,  is  the  effect  of 
convection  in  the  melt  on  thermal  gradients  and  composition  in  the  solid  phase 
(for  example,  see  Derby  [2]).  By  accurately  representing  the  melt  interface, 
fluid  dynamics,  and  heat  transfer,  one  can  better  predict  crystal  characteristics 
based  on  the  furnace  geometry  and  thermal  environment.  In  addition,  because 
the  growth  process  is  lengthy,  there  is  the  possibility  of  active  control  of  certain 
crucible  factors,  such  as  applied  magnetic  fields  and  time-dependent  heating. 

Mr.  Johansen's  approach  is  based  on  previous  work  with  “Cartesian  grid” 
and  “volume-of-fluid”  methods  for  unsteady  fluid  dynamics.  These  methods 
represent  the  problem  on  a  square  finite  difference  grid  and  apply  standard  dis¬ 
cretization  techniques  to  the  interior  of  the  domain.  Near  irregular  boundaries, 
on  a  smaller  number  of  points,  special  algorithms  are  used  to  advance  the  solu¬ 
tion,  taking  the  local  geometry  into  account.  These  Cartesian  grid  methods  were 
recently  applied  to  inviscid  compressible  flows  by  Pember,  et  al.[3].  Time  step 
constraints  (due  to  small  cells)  were  overcome  by  redistributing  excess  fluxes  to 
larger  neighboring  cells,  while  also  maintaining  conservation.  The  approach 
also  includes  AMR,  local  Adaptive  Mesh  Refinement  of  the  finite-difference 
grid,  based  on  work  done  by  Berger  and  Colella  [4];  results  for  two-  and  three- 
dimensional  regions  compare  well  with  other,  more  expensive  computations.  A 
similar  approach  has  been  developed  for  incompressible  flows  by  Almgren  et 
al.  [13],  which  included  efficient  solution  of  a  Poisson  equation,  with  no-flux 
boundary  conditions,  on  the  irregular  domain.  Volume-of-fluid  algorithms  have 
been  applied  to  a  number  of  problems  in  compressible  flow.  Chern  and  Colella 
[5]  implemented  conservative  interface  tracking  for  compressible  flow,  using  a 
Simple  Line  Interface  Calculation  (SLIC),  to  represent  the  boundary  between 
different  gases.  When  combined  with  AMR,  the  method  produced  excellent  cor¬ 
relation  with  experiments  of  shocks  impinging  on  gas  interfaces  [6],  even  with 
the  first-order  representation  of  the  front.  Since  then,  similar  volume-of-fluid 
algorithms  have  been  developed  for  multifluid  compressible  flows  [7],  even  in¬ 
cluding  combustion  [8],  and  complicated  physical  domains  [9]. 

These  papers  demonstrate  that  modeling  the  melt  motion  is  certainly  within 
the  capabilities  of  a  volume-of-fluid  method.  The  largest  hurdle  in  applying  it 
to  crystal  growth  is  developing  an  algorithm  for  tracking  the  interface,  while 
maintaining  conservation  of  energy.  Once  this  is  accomplished,  a  multi-step, 
predictor-corrector  algorithm  can  be  used  to  advance  the  solution  in  time,  in 
which  both  the  melt  interface  and  container  geometry  will  be  represented  using  a 
hybrid  front-tracking  method  derived  from  work  mentioned  above.  Combining 
the  Stefan  problem  with  an  integrator  for  the  Boussinesq  equations  for  the 
melt  (  similar  to  recent  progress  for  incompressible  flows  [11,  13,  14]  ),  will  be 
the  last  step  for  the  base  algorithm.  Adding  an  adaptive,  finite-difference  grid 
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hierarchy  will  provide  local  refinement  of  interfaces,  and  help  verify  the  method. 
Implementation  of  the  complex  algorithm  is  simplified  with  the  use  of  hybrid 
and  Fortran  libraries  developed  at  Lawrence  Livermore  National  Labs. 

The  details  of  the  method  are  based  on  the  papers  mentioned  above,  with 
adaptations  and  improvements  for  the  Stefan  problem.  Our  front  reconstruc¬ 
tion  algorithm  will  be  based  on  work  done  by  Pilliod  and  Puckett  [10].  Their 
approach  provides  a  method  of  approximating  an  interface,  based  on  a  linear, 
least-squares  fit  to  the  local  array  of  volume  fractions.  The  algorithm  is  shown 
to  be  second-order  accurate  in  space,  even  when  applied  to  complicated  shapes. 
For  the  Stefan  problem,  the  interface  moves  normal  to  itself,  with  a  speed  pro¬ 
portional  to  the  local  temperature  gradient.  Chern  and  Colella  [5]  implemented 
shock  tracking  for  compressible  flow,  which  is  directly  applicable  to  the  motion 
of  the  front  in  the  Stefan  problem.  By  coupling  the  reconstruction  algorithm 
with  an  operator-split  version  of  this  front-tracking  scheme  (  modified  to  be 
minimum  and  maximum  preserving  ),  we  will  be  able  to  advance  the  phase 
boundary  in  time.  Conservation  of  energy  is  enforced  by  adding  equivalent  heat 
sources,  due  to  enthalpy  of  formation,  where  necessary. 

Solving  the  temperature  equation  on  the  irregular,  changing  domain  will  be 
the  next  step.  First,  the  half-time  velocities  will  be  found  using  an  algorithm 
similar  to  those  in  [13].  An  intermediate  MAC  projection  guarantees  that  the 
edge- centered  velocities  have  zero  discrete  divergence,  making  them  ideal  for 
computing  a  time- centered,  temperature  convection  term.  This  will  then  be 
used  in  the  right-hand  side  of  the  temperature  equation,  along  with  the  en¬ 
thalpy  of  formation  mentioned  earlier.  The  melting  point  of  the  substance  will 
require  a  Dirichlet  boundary  condition  for  temperature  on  the  interface;  how¬ 
ever,  prescribing  the  front  motion  is  equivalent  to  specifying  the  temperature 
gradient  at  the  interface,  too.  This  conflict  will  be  surmounted  by  ignoring  the 
constraint  on  the  temperature  gradient;  it  is  determined  by  fitting  a  quadratic 
through  values  interpolated  from  full  cells  near  the  interface,  using  the  normal 
and  midpoint  of  the  reconstructed  front.  The  implicit  discretization  of  the  diffu¬ 
sion  equation  will  be  solved  using  multi-grid  iterations  to  accelerate  convergence 
[11]  [12]  .  Once  the  new  temperature  field  is  determined,  the  energy  loss  due 
to  a  mismatch  in  gradients  will  be  added  into  the  temperature  equation  on  the 
next  time  step. 

Finally,  the  front  position  and  temperature  fields  will  be  used  in  the  inte¬ 
gration  of  the  Boussinesq  equations  for  flow  in  the  melt.  The  remainder  of  the 
algorithm  is  similar  to  that  presented  in  [13]  and  [14].  The  main  differences 
will  be  the  inclusion  of  viscosity,  which  demands  no-slip  boundary  conditions  at 
all  interfaces.  The  viscous  forces  will  be  treated  similarly  to  the  heat  equation, 
with  the  advection  terms  and  body  forces  (due  to  gravity  and  externally  applied 
magnetic  fields)  serving  as  sources  in  the  Crank-Nicholson-type  discretization. 
In  the  final  step,  the  velocity  field  will  be  projected  onto  its  divergence-free  part. 

Progress  to  date  has  focussed  on  accurate  integration  of  the  Stefan  prob¬ 
lem,  and  the  heat  equation  with  Dirichlet- type  boundary  conditions.  First, 
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Mr.  Johansen  has  thoroughly  tested  a  one-dimensional  method  for  the  classical 
Stefan  problem,  which  represents  a  balance  between  heat  conduction  and  en¬ 
thalpy  of  formation.  The  interface  gradient  is  calculated  using  a  quadratic  fit 
to  nearby  temperature  values;  after  the  front  is  advanced,  conservation  is  en¬ 
forced  by  adding  appropriate  source  terms.  The  quadratic  temperature  fit  and 
lagged  conservation  have  been  shown  to  be  stable  for  certain  discretizations; 
results  demonstrate  second-order  convergence  in  space,  and  first-order  accuracy 
in  time,  regardless  of  Stefan  number.  Second,  the  solution  of  Dirichlet-type 
Poisson  problems  on  irregular,  two-dimensional  domains,  provides  part  of  the 
time-stepping  algorithm  for  the  Stefan  problem  in  more  than  one  dimension. 
The  calculation  of  interface  gradients,  bcised  on  quadratic  fits  to  interior  val¬ 
ues,  has  been  demonstrated  as  globally  second-order  accurate,  even  for  front 
representations  containing  arbitrarily  small  cells. 

Mr.  Johansen  will  have  a  working  code  by  December,  1996.  Adaptive  so¬ 
lution  of  the  Stefan  problem  in  multiple  dimensions  will  be  finished  by  June, 
1996;  it  will  provide  the  crucial  step  needed  to  incorporate  the  front  motion  into 
a  Cartesian  grid,  incompressible  Navier-Stokes  solver,  which  is  being  developed 
at  both  UC  Berkeley  and  LLNL,  and  will  be  completed  by  September,  1996. 
The  remaining  three  months  will  be  spent  testing  the  algorithm,  by  verifying 
it  against  experimental  and  computational  results;  externally-applied  magnetic 
fields  might  be  included  in  the  model  in  Spring  of  1997.  This  work  will  be  a 
major  step  in  applying  finite-difference  algorithms  to  real  engineering  problems. 
The  simulation  of  semiconductor  substrate  crystal  growth  will  be  a  new  appli¬ 
cation  of  the  Cartesian  grid  method,  with  front  tracking,  convection  of  heat  and 
melt  comi>osition,  and  adaptive  refinement  of  the  finite-difference  grid  around 
significant  features. 
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